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A model for epidemics on an adaptive network is considered. Nodes follow an SIRS (susceptible- 
infective-recovered-susceptible) pattern. Connections are rewired to break links from non-infected 
nodes to infected nodes and are reformed to connect to other non-infected nodes, as the nodes that 
are not infected try to avoid the infection. Monte Carlo simulation and numerical solution of a mean 
field model are employed. The introduction of rewiring affects both the network structure and the 
epidemic dynamics. Degree distributions are altered, and the average distance from a node to the 
nearest infective increases. The rewiring leads to regions of bistability where either an endemic or 
a disease- free steady state can exist. Fluctuations around the endemic state and the lifetime of the 
endemic state are considered. The fluctuations are found to exhibit power law behavior. 



I. INTRODUCTION 

The study of recurrent epidemics has a long history 
and many models, both deterministic and stochastic, 
have been considered. Deterministic models have been 
used since the time of Bernoulli and have explained some 
of the mechanisms in the spread of the infectious dis- 
eases. However, deterministic models are not sufficient 
to account for some of the important stochastic dynam- 
ics, such as extinction [2, 3] and sustained fluctuations 
4|. From the general theory of finite Markov chains 
5|, it was shown that in stochastic models the proba- 
bility of extinction is equal to one in the asymptotic time 
limit. Numerical @, U [1] and analytic ^Qj comparisons 
of stochastic and deterministic models have been per- 
formed. The numerical results hold for very small ampli- 
tude noise as well as real finite noise. Deterministic SIS 
or SIRS models result in an equilibrium endemic presence 
of infectives for an appropriate choice of parameters. It is 
clear that stochastic effects may result in very different 
dynamics from deterministic models, particularly when 
fluctuations and/or extinction occur. 

More recently, the study of fluctuations and the spread 
of simple models of epidemics have been simulated on 
large networks [13, [HI El, [H, • In almost all of these 
network models, the epidemic propagates on a fixed net- 
work. The epidemic dynamics is typically studied as an 
SIS or SIR model, in which the population is large and 
isolated. In addition to the dynamics on such fixed ar- 
chitectures, controls based on vaccination have been con- 
sidered as well [H, [31 ■ Several recent models have con- 
sidered epidemics on a network that changes structure 
dynamically according to rules that do not depend on 
the nodes' epidemic status [TtI. [isj. 

In contrast to the models of a static network or models 
with externally applied changes in structure, a new class 
of models based on endemic SIS populations on an adap- 
tive network has been recently introduced . Changes 
to the network structure are made in response to the epi- 



demic spread and in turn affect future spreading of the 
epidemic. Here, the new parameter is one that describes 
the rewiring rate of the network, which is controlled by 
the fraction of susceptible (S)-infective (I) links. The net- 
work alters dynamically when there are contacts between 
S and I, and social pressures (the desire to avoid illness) 
rewire the contacts to be instead between S and S. Infec- 
tions are reduced due to isolation, and a new phenomenon 
occurs: for appropriate choices of parameters, bistability 
between the disease free equilibrium and endemic state 
has been observed. This is in contrast to static networks 
in a large population, where there is typically only a sin- 
gle attracting endemic or disease free state. 

In this paper, we introduce a recovered, immune class 
and consider this slight generalization of the SIS model 
on an adaptive network. We examine the structure of 
the network and the dynamics of the fluctuations of the 
epidemic. Our approach is to combine Monte Carlo sim- 
ulations and stochastic mean held models for epidemic 
evolution on evolving networks. The layout of the paper 
is as follows. We introduce the model in Section [H] and 
present its bifurcation structure in Section IIIII Proper- 
ties of the network structure are discussed in Section [TVl 
We discuss dynamical properties of the system, including 
fluctuations and lifetimes of the states, in Section IVl 



II. MODEL 

We study an SIRS (susceptible-infective-recovered- 
susceptible) model on an adaptive network. Epidemic 
dynamics on the nodes is as follows. The rate for a sus- 
ceptible node to become infected is piV/^nbn where iV/^nbr 
is the number of infected neighbors the node has. The 
recovery rate for an infected node is r. A recovered node 
becomes susceptible again with rate g, which we deflne 
as the re-susceptibility rate. 

While the epidemic spreads, the network is also be- 
ing rewired adaptively. If a link connects a non-infected 
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node to an infected node, that link is rewired with rate 
w to connect the non-infected node to another randomly 
selected non-infected node. Self links and multiple links 
between nodes are disallowed. 

In examining steady state solutions, it is sufficient to 
fix one of the rates, as time may be rescaled accordingly. 
For this reason, we fix r = 0.002 throughout this paper. 

We performed Monte Carlo simulations of this model 
on a system with N — lO"' nodes and K = 10^ links. 
(Larger system sizes with the same node to link ratio 
were also considered. The major results of this paper 
do not depend strongly on system size.) In each Monte 
Carlo step (MCS), we randomly select N nodes and M 
links, where M is the number of links that may poten- 
tially rewire (susceptible-infected and recovered-infected 
links), and the links are selected from the pool of links 
that may rewire. Initial conditions are constructed in one 
of two ways. We either generate a random (Erdos-Renyi) 
graph of susceptibles and convert a fraction / of them to 
infectives, or we use the final state of a previous run as 
an initial condition. Transients are discarded and simu- 
lations run long enough that the initial conditions do not 
affect the results. 

Following [l^, we also developed a corresponding 
mean field model for the system. The mean field model 
tracks the dynamics of both nodes and links. Pa denotes 
the probability of a node to be in state A, where A is ei- 
ther S (susceptible), I (infected), or R (recovered). Pab 
denotes the probability that a randomly selected link con- 
nects a node in state A to a node in state B. We obtain 
the following mean field equations for the evolution of 
the nodes: 
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For example, in the first equation, recovereds are con- 
verted to susceptibles with rate q, and infection spreads 
with rate p along each susceptible-infected link. Rewiring 
does not appear directly in the node equations, since 
rewiring operates on links, but it affects the system im- 
plicitly through the number of susceptible-infected links 
[KPsi). We next write a system of mean field equa- 
tions for the links. To close the system, we follow [1^ 
and make the assumption for three point terms that 
Pabc ~ PabPbc I Pb- This assumption leads to the 
following system of equations for links: 
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The ordinary differential equations of the mean field 
model can be integrated easily with any well known nu- 
merical integration technique. We choose initial condi- 
tions so that we are near an endemic state. We also note 
that the model does support solutions with negative val- 
ues, but these are unphysical, and so we ignore them. In 
the case of stochastic simulations, we have considered the 
effects of both internal fluctuations, modeled as multi- 
plicative noise, as well as external fluctuations, modeled 
as additive noise. We use a fourth order Runge-Kutta 
solver for each of these cases to generate stochastic stim- 
ulations of the mean field. We also tracked the steady 
states as a function of parameters using a continuation 
package [20| . 



III. BIFURCATION STRUCTURE 

We first consider the steady state bifurcation struc- 
ture of the model. In Figure [Tl we show examples of 
the average infected fraction versus the transmission rate 
p. Two steady states can occur, a disease-free state and 
an endemic state. In the absence of rewiring (Fig. [T^), 
the disease-free state loses stability for very small trans- 
mission rates, and only the endemic state is observed at 
larger p values. When rewiring is introduced (Fig. [T|d), 
the disease-free state is stabilized for larger p values. A 
region of bistability, in which both endemic and disease- 
free states are observed, now occurs. Bistability was ob- 
served previously for the SIS model [19]. Monte Carlo 
simulation points in Fig. [1] were computed as follows. 
To locate the upper branches (endemic state), we swept 
p from larger to smaller values, using the final state of 
each run as the initial state for the next run. 5 x lO'^ 
MCS of transients were discarded, and the steady state 
was averaged over 10^ MCS. To locate the lower branch 
(disease free state) , we began with a randomly connected 
network in which a fraction / of the nodes were infected, 
while the rest were susceptible. The system was simu- 
lated for 2 X 10" MCS. / values between 0.025 and 0.9 
were tried, and at least 5 runs were done for each p value. 
If the epidemic died out in any of the runs, the disease 
free state was considered stable. Due to the stochastic 
nature of the Monte Carlo simulation, and because the 
disease-free state is absorbing, stability designations are 
uncertain. It is difficult to distinguish a weakly unstable 
state from a weakly stable state with a short lifetime. 
Lifetimes of the endemic state are considered in more 
detail in Section IVl 

The results in Figure [T] show fairly good agreement be- 
tween the mean field approximation and the Monte Carlo 



3 



0.008 




0.008 



FIG. 1: Average infected fraction vs. transmission rate p. 
(a) Static network, w = Q. (b) Rewired network, w — 0.04. 
Black dots: Monte Carlo simulations; solid gray lines: mean 
field solution (stable branches); dashed gray lines: mean field 
solution (unstable branches), q = 0.0016, r = 0.002. 



simulation of the full system, and we typically see this 
level of agreement in the steady state values, although 
as we discuss later, the stability and type of bifurcation 
sometimes differs. Using the mean field model, we next 
explore the bifurcation structure of the system for a wider 
range of parameters. 

An interesting property of the steady state instabili- 
ties appears when one considers each of the steady state 
bifurcation points of the mean field equations. (There 
do exist branches of periodic orbits, but since they oc- 
cur within a very small range of parameters, we ignore 
them in this paper. They will be treated elsewhere.) If 
the re-susceptibility rate q is held fixed and a bifurcation 
diagram constructed, we find the existence of at least 
two distinct regimes for different q values, illustrated in 
Figs. and [31 The instabilities appear as a transcriti- 
cal bifurcation from the disease free state steady state, 
a saddle-node bifurcation of endemic steady states, and 
a Hopf bifurcation, from which a branch of subcritical 
unstable periodic orbits emanates (not shown). 

In the case where q = 0.0064, depicted in Fig. [21d for 
w — 0.04, we show a typical bifurcation plot of stable and 
unstable steady states as the infection rate p is increased. 
The dashed lines represent the unstable steady states. 
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FIG. 2: (a) Two parameter plot for w and p of the bifurcation 
points for steady states when q — 0.0064. The heavy dashed 
line is the line of saddle-node points (limit points). The solid 
line denotes the Hopf bifurcation points, and light dashed 
line denotes the transcritical bifurcation points (BP) . (b) A a 
bifurcation diagram of the infective fraction as a function of 
p, where w = 0.04. The squares denote the saddle-node point 
and transcritical point. Dashed lines are unstable branches. 



while the solid lines depict the attractors. For low values 
of p the disease free steady state is stable. Tracking along 
the disease free branch, at a critical value of p, an unsta- 
ble branch (subcritical) of endemic steady states appears. 
The endemic branch becomes stable at the saddle-node 
point. There exists a clear region of bistability with co- 
existing endemic and disease free states for a range of p. 
If we now vary the parameter w and track each bifurca- 
tion curve, we obtain the result in Fig. [5^. We describe 
the bifurcation regions in detail for w larger than 0.1. In 
region I, we have only a stable disease free equilibrium. 
As we cross into region II for large w, the disease free 
equilibrium is stable, and there exists an unstable en- 
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demic state. Region III exhibits bistability between the 
disease free equilibrium and endemic state, and region 
IV has just a stable endemic equilibrium. We note that 
at w — 0.04, we have the simple saddle-node transition 
depicted in Fig. [2Jd, since there is no Hopf bifurcation to 
periodic cycles at that particular w value. 

The previous discussion presented a case for re- 
susceptibility rate q where the limit (saddle-node) and 
Hopf bifurcation branches are close to each other. The 
distinction between the saddle-node and Hopf branches 
can be seen more easily if q is lowered to 0.0016, as 
shown in Fig. [3^. In Fig. [3^, for sufficiently large w, 
as p is increased in region I the system first undergoes 
a limit point bifurcation, and then a Hopf bifurcation as 
it passes through region II. The Hopf curve is actually a 
closed isola in two parameters. The limit point here is 
a saddle-saddle point, where a steady state having a two 
dimensional unstable manifold connects to a steady state 
with a one dimensional unstable manifold. In both cases, 
we have bistable behavior for w sufficiently large, but the 
region of bistability is much smaller since the Hopf and 
transcritical branches are closer together for this value of 
q (region HI). 

For w = 0.04, the mean field endemic steady state loses 
stability in a saddle-node bifurcation for q = 0.0064 and 
in a Hopf bifurcation for q = 0.0016. In our discussion 
of fluctuations in the endemic state in Section |Vl we will 
refer to q — 0.0064 because the saddle-node bifurcation 
structure best corresponds to the scaling of fluctuations 
that we observe in Monte Carlo simulations of the full 
system. 
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IV. NETWORK GEOMETRY 
A. Degree distributions 

Rewiring leads to significant alterations in the network 
structure. We first consider the degree distribution. Fig- 
ure m shows degree distributions for each type of node 
in the absence (jl^) and presence ([Us) of rewiring. In 
Fig. [4)3, we averaged over 3 x 10"* MCS. In Fig. [li, the 
network is static, so we averaged over 10 separate runs 
to obtain better statistics. 

Using mean field ideas, we can understand the recov- 
ered and susceptible degree distributions as nodes flow 
from infected to recovered to susceptible. We outline 
the calculation briefly here. However, this approach does 
not accurately predict the degree distribution for infected 
nodes because correlations play a more important role for 
these nodes, as we will explain below. For this reason, we 
cannot write a self-consistent set of equations for the de- 
gree distributions that could be solved without inputting 
simulation results. 

Let dx,n be the number of nodes in state X (either 
S, I, or R) with degree n. Recovered nodes originate 
when infectives recover and are lost when they become 
susectible again. The degree of a recovered node can only 



FIG. 3: (a) Two parameter plot of the steady state bifurca- 
tion points when q — 0.0016. The heavy dashed line is the line 
of limit points. The solid line denotes the Hopf bifurcation 
points, and the light dashed line denotes the transcritical bi- 
furcations, (b) A bifurcation diagram of the infective fraction 
as a function of p, where w — 0.04 (same bifurcation diagram 
in Fig.[lJ. 



increase, as other susceptible and recovered nodes wire to 
connect to it. This leads to the following equations for 

dR.n- 



dt 
d 

di 



(dRfi) = rdifi - qdn. 
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(10) 



where k is the average rate for nodes to rewire to a given 
non-infected node, k is given by the ratio of the total 
rewiring rate to the number of potential target nodes: 



wKjPsi + Pir) 
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(11) 



Given the degree distribution of the infectives d/^„ and 
the probabilities appearing in k, Eqs. [TUI can be solved 
for the degree distribution of the recovereds. 
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FIG. 4: Degree distributions from Monte Carlo simulation for 
p = 0.002, q = 0.0016, r = 0.002. (a) Static network, w = 
0. All node types have Poissonian degree distribution, (b) 



Rewired network, w = 0.04. Solid gray: 
recovereds; black: susceptibles. 



infectives; dashed: 



Susceptible nodes originate when recovereds become 
susceptible again, and they may be lost when they be- 
come infected by a neighbor. As with the recovereds, the 
degree of a susceptible increases due to rewiring. Thus 
the time evolution of the degree distribution for the sus- 
ceptibles can be written as 



dt 
d_ 

di 



(ds.o) = qdR,Q - kds.o 

ids,n) = qdR,n - k'nds,n - kds,n 
+kds,n-i for n > 



(12) 



where k' is the infection rate per link into a suscepti- 
ble node. We assume that k' is independent of degree, 
which we know from simulations is approximately correct 
(cf. Fig.O, and write 



fc' =p 



SI 



Psi + PsR + 2Pss 



(13) 



where the fraction is the ratio of the number of links that 
can transmit infection to the total number of links into a 
susceptible. As with the recovereds, the steady state de- 
gree distribution for susceptibles can be computed from 
Eq. [121 The predicted degree distributions for suscepti- 
bles and recovereds are overlaid on the actual distribu- 
tions in Fig. 0)3, using Monte Carlo simulation averages 
for the infective degree distribution and the node and 
link probabilities. (Note that the node and link prob- 
abilities could instead be obtained from the mean field 
system.) Deviations between the prediction and simula- 
tion are smaller than the width of the curves in Fig. [3]d, 
so they are indistinguishable. 

The degree of infected nodes, however, cannot be pre- 
dicted by this approximate procedure. We might expect 
that 



— k"dj^i — rdjo 

= k"{n + l)di,„+i-k"ndi,n (14) 
+k'nds,n — rdi^n for n > 
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is the per link rewiring rate for links connecting to an in- 
fective (i.e., the ratio between links that can potentially 
rewire and total links into infectives). Figure [5^ com- 
pares the actual degree distribution for infectives with 
the distribution predicted by Eq. 1141 The number of low 
degree infectives is significantly over-predicted. This oc- 
curs because the mean field approximation in Eq. [15] is 
not accurate for infectives. Figure [5)3 shows the frac- 
tion of infected neighbors that a node has, depending 
on its degree and disease status. Results are averaged 
over 3 x 10^ MCS. Low degree infected nodes tend to 
have a much higher fraction of neighbors that are also in- 
fected, due to transmission of the disease. Once infected, 
these neighbors will not rewire away until recovered, so 
the rewiring rate per link is smaller than one would ex- 
pect from the mean field k" in Eq. [T5| Predicting the 
infective degree distribution accurately would require a 
theory that accounts for these correlations in infection 
status of neighboring nodes, which is beyond the scope 
of the present work. 



B. Distance from an infective 

We next consider the distribution of distances from a 
given node to the nearest infective. These distances are 
of interest because they relate to the number of hops the 
disease must make in order to reach an uninfected indi- 
vidual. The disease cannot propagate through recovered 
nodes until they become susceptible again, so the dis- 
tance from the nearest infective does not necessarily cor- 
respond to a path for disease propagation. However, we 
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FIG. 5: (a) Actual (black) and predicted (dashed gray) degree 
distributions for infecteds. (b) Average fraction of neighbors 
that are infected vs. degree. Solid gray: infecteds; dashed: 
recovereds; black: susceptibles. p = 0.002, q — 0.0016, r — 
0.002, w = 0.04. 



note that rewiring acts only on links to infectives, and 
thus the chains of susceptible and recovered links that 
this metric identifies will persist until the infection leads 
to their interruption. Figure [B^ shows the distribution of 
distances from the nearest infective in the presence and 
absence of rewiring. To display the effect on the network 
geometry alone, rather than on the steady state number 
of infections as well, we have used a smaller q value for 
the w = case so that the total number of infectives is 
approximately the same in both curves. Results were av- 
eraged over 3 x 10"* MCS, sampled every 100 MCS after 
removing transients. Rewiring significantly decreases the 
number of nodes that are directly connected to an infec- 
tive. However, despite the rewiring, only a small fraction 
of nodes are fully disconnected from the infection. 

Figure [6b shows a semilog plot of the w = 0.04 case 
to display the tail of the curve at larger distances. An 
approximation based on random networks is also shown. 




distance 



FIG. 6: Distribution of distances from the nearest infected 
node, oo indicates nodes that are completely disconnected 
from an infective, (a) Solid black: with rewiring {w — 0.04, 
p = 0.002, q — 0.0016); dashed gray: no rewiring {w — 0, 
p = 0.002, q — 0.009). (b) Black: rewired case (same as in 
(a)); gray: distribution for random graphs, as described in 
text. 



Given how far from Poisson the degree distributions are 
when the network is rewired, it is somewhat surprising 
that the form of the decay in the distribution of distances 
can be predicted from random networks. Beginning with 
S, I, and R nodes in numbers matching that observed 
in the average of Monte Carlo simulations, we generated 
1000 random networks and added randomly selected links 
until the number of S-S, S-I, S-R, etc., links also matched 
that from the average of Monte Carlo simulations. As 
Figure [61d shows, the distribution of distances for these 
random networks decays in the same way as it does in an 
adaptive network. Thus the form of the distribution of 
distances depends mainly on local dynamics (node and 
link dynamics) rather than on the details of longer range 
correlations. The main difference is that the adaptive 
network has some nodes fully disconnected from infected 
components, while the random networks do not. Most of 
these nodes are recovereds with degree 0, which appear 
when infectives of degree recover. Each type of node in 
the random networks has a Poisson degree distribution, 
so they do not generally have nodes of degree 0. 
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V. FLUCTUATIONS AND OTHER DYNAMICS 

In the previous sections, we have considered steady 
states and long time averages of network properties. We 
next consider fluctuations and dynamical properties of 
the endemic state. 



A. Fluctuations near bifurcation point 

Near the bifurcation point where the endemic state 
loses stability, the number of infectives has larger fluc- 
tuations due to noise overcoming weak attracting forces 
[2H . Fluctuations in the SIRS model are significantly 
larger than those in the previously studied SIS model. 
We quantify the fluctuations by computing the standard 
deviation divided by the mean for long time series in both 
the Monte Carlo and mean field simulations. In Fig. [3 
we plot the fluctuations as the infection rate p is swept 
towards the bifurcation point. Monte Carlo results were 
computed from 5 x 10^ MCS time series sampled every 
10 MCS, except for the two smallest p values, for which 
shorter time series were used due to the shorter lifetimes 
of these states. All time series were longer than 10^ MCS. 

For comparison, the mean field equations can also be 
considered near equilibrium in stochastic form. In gen- 
eral, near equilibrium fluctuations can be modeled as 
additive noise [l^, and we do so here. (We do note 
that multiplicative noise effects generate similar results 
to those reported for additive noise.) We assume the 
mean field is of the following form: 



X' = F(X)+er;(t), 



(16) 



where F(X) is the mean field system in Eqs. [T][51 and 
(?7(t)77(t')) = S{t — t'). e is the noise strength, or ampli- 
tude. We have considered both additive noise and multi- 
plicative noise cases in the simulations of the stochastic 
attractors near the endemic state and have computed 
the standard deviation divided by the mean as described 
above for 10 random initial conditions near equilibrium 
and 10 realizations. 

We recall from Figs. [2] and [3] that depending on the 
value of the re-susceptibility rate q, the bistability regions 
III are vastly different. Specifically, for q = 0.0016, we 
saw that for sufficiently large wiring rates, the saddle- 
saddle and Hopf bifurcation branches were well sepa- 
rated, whereas for q — 0.0064, the branches were very 
close for small values of re- wiring rate w. In Monte Carlo 
simulation, we have not observed the Hopf bifurcations or 
stable periodic oscillations seen in the mean field, even for 
system sizes as large as 4 x 10^ nodes. Although the value 
oi q — 0.0064 we use in the mean field fluctuation study 
is different from that used in the Monte Carlo, the local 
bifurcation structure is similar when w = 0.04, in that 
it is a true saddle-node bifurcation point. This has been 
checked by examining the local linear vector field at the 
saddle-node point in question. Therefore, although the 



mean field has different q value, the bifurcation structure 
is equivalent to that observed in Monte Carlo simulation, 
so we use q — 0.0064 in the fluctuation study. 

The computations reveal that the fluctuations exhibit 
power law scaling, as shown in the log- log plots of Fig.[7|3 
and[7ti. On the horizontal axis, we plot ln(j3 — pc), where 
Pc is the critical point at which the endemic state loses 
stability. For the mean field, the bifurcation point is 
known exactly by examining the eigenvalues of linearized 
vector field at the steady state. We approximate the 
Monte Carlo bifurcation point Pc as the value that pro- 
duces the most linear plot. Although both cases have 
power law scaling, the exponents are different: —0.59 for 
Monte Carlo and —0.27 for mean field. The scaling expo- 
nent for the full system depends on the number of nodes. 
Whether it will approach the mean field value in the limit 
of infinite system size is a subject for further study. We 
note that fluctuations near a Hopf bifurcation point, as 
occurs in the mean field for q — 0.0016, would produce a 
very different form of scaling from the saddle-node case. 
One of the reasons for the difference is that the instabil- 
ity is then two dimensional and underdamped [23j . This 
is known to cause very different scaling laws in generic 
problems, which can be much slower [24j |. 

To motivate the power law scaling of the fluctuations, 
we consider scaling near a generic saddle-node bifurca- 
tion. The simplest generic case of a saddle-node bifur- 
cation for equilibrium points comes from solving for ze- 
roes of the vector field at a parameter value where one 
eigenvalue of the Jacobian passes through zero. Standard 
normal form analysis allows one to consider the generic 
problem of a saddle-node bifurcation. In one dimension, 
the stochastic differential equation of a saddle-node bi- 
furcation may be modeled as 



dxt = (a — x1)dt + a * dWt- 



(17) 



The parameter is a, and we suppose noise is additive. 
Since noise in general may cause a "shift" in parame- 
ter values where the saddle-node point disappears, we 
assume that the noise near the bifucation is sufhciently 
small, where dW/dt is a white noise term, and dW is a 
Brownian increment. 

We further assume that we are always near the attract- 
ing branch of the saddle-node, so we are in a near equi- 
librium setting. Such an assumption allows us to exam- 
ine the stationary probability density function (PDF) of 
the stochastic dynamics by employing the Fokker-Planck 
equation near steady state. For the stochastic differential 
equation, Eq.[T71 the PDF is well known 2l| and is given 

by 



p(a, X, cr) = iVe 



- /\r^=2(aa;-a;V3)/'T2 



(18) 



Here is a normalization constant. We compute the 
first and second order moments directly using Eg . [T51 and 
then take the ratio of the standard deviation to the mean. 
Since a = is the value of the saddle-node point, we ex- 
amine the fluctuations in the neighborhood of that value. 
The results, shown in Fig. [51 display power law scaling. 
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FIG. 7: Fluctuations in infectives (standard deviation divided by mean) vs. infection rate p near the bifurcation point: (a) 
Monte Carlo, (b) mean field. Curves are to guide the eye. Log-log plots (data points with best fit lines) show power law scaling 
for both Monte Carlo (b) and mean field (d). Monte Carlo parameters: q = 0.0016, w — 0.04, r — 0.002. Mean field parameters: 
q = 0.0064, w = 0.04, r = 0.002, e = 0.0001. 



B. Delayed outbreaks 

We next consider phase relationships between the fluc- 
tuating variables. We tracked the number of infectives 
in the system at each time point as well as the number 
of nodes that neighbored an infective (i.e., the number 
of SI and IR links). In the rewired system, fluctuations 
in the number of infectives lagged behind fluctuations in 
the number of infective neighbors that arc not themselves 
infected, as shown in Fig. [5^. 

Both mean field and Monte Carlo simulations of the 
full system displayed this effect, and we studied its de- 
pendence on the rewiring rate. Monte Carlo simulations 
were sampled every 1 MCS for 3 x 10'' MCS after discard- 
ing transients. The mean field was sampled every 1 time 
unit for 5 x 10^ units. Additive noise was included in 
the mean field equations with noise strength e = 0.0001. 



Cross correlations between the infectives and the infec- 
tive neighbors were computed for varying shifts between 
the time series, and the lag maximizing the cross corre- 
lation was identified. As shown in Fig.[5]3, rewiring leads 
to increasing lag times and delayed outbreaks. We also 
computed time lags for the mean field model with multi- 
plicative noise and found the same trend of increasingly 
delayed outbreaks with larger rewiring. 



C. Lifetime of endemic steady state 

The final dynamic effect we consider is the lifetime of 
the endemic steady state. Because the system is stochas- 
tic and the disease free state is absorbing, all parameter 
values will lead to eventual die out of the disease in the in- 
finite time limit. These lifetimes become shorter and die 
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FIG. 8: Fluctuation size of a generic saddle-node bifurcation 
as a function of bifuration parameter a near the bifurcation 
point using Eq. 1181 The noise level used is 0.005. 



FIG. 10: Dependence of endemic state average lifetimes on 
infection rate p. Points: Monte Carlo simulations; line: best 
fit line, q = 0.0016, r = 0.002, w = 0.04. See text for details. 
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FIG. 9: Delayed outbreaks due to rewiring, (a) Monte Carlo 
time series. Black: infectives; gray: neighbors of infectives. 
Curves are scaled in arbitrary units for comparison of peak 
times, p = 0.0065, w = 0.09, q = 0.0016, r = 0.002. (b) Time 
in MCS by which infectives lag behind infective neighbors 
vs. rewiring rate. Solid black: Monte Carlo; dashed gray: 
mean field, p = 0.0065, q = 0.0016, r = 0.002, mean field 
noise strength e = 0.0001. 



out is more easily observed in the bistable regime, near 
the bifurcation point where the endemic state has weak 
stability. We measured the dependence of the lifetimes 
on the infection rate p. For each p value, we prepared 
a steady state initial condition and computed multiple 
duplicate runs to obtain a distribution of lifetimes. (We 
computed 100 duplicate runs for all but the two high- 
est p values and over 25 runs for the two highest.) We 
then calculated the average lifetime T for each p. For 
a generic saddle-node bifurcation in one dimension, it 
is expected that In T is a linear function of [p — Pof^"^ ^ 
where po is the bifurcation point [H, [1^ . We obtained 
spurious results for our system with 10* nodes, possibly 
because the small system size led to rapid die out and we 
were unable to obtain good statistics near the bifurcation 
point. Away from the bifurcation point, the system has a 
weakly damped oscillatory component and behaves like 
a focus. By switching to a system with 4 x lO'* nodes 
and 4 x 10^ links, we were able to run longer simulations 
closer to the bifurcation point and operate in a regime 
where only one dimension mattered and the oscillations 
could be ignored. Preliminary scaling results are shown 
in Figure [TOl We used the bifurcation point po that gave 
the best fit line, which led to an R value of 0.99. The 
scaling results appear consistent with expectations, but 
more study and better statistics are needed. 

VI. CONCLUSIONS AND DISCUSSION 

We have explored the stable states, network proper- 
ties, and dynamics of an SIRS model on a network with 
adaptive rewiring. As with the SIS model studied previ- 
ously by Gross et al. [l^ , the rewiring leads to bistability 
of the endemic and disease free states. A mean field ver- 
sion of the model predicted the steady states with good 
numerical accuracy and was also valuable in studying the 
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fluctuations of the system, with the caveat that one must 
be near the appropriate type of bifurcation in the mean 
field to obtain corresponding results. With the addition 
of the recovered class and re-susceptibility rate, we can 
control the width of the bistability region by manipulat- 
ing the location of the bifurcation points. 

The fluctuations in the infectives near the bifurca- 
tion point showed power law scaling. This agreed with 
mean field results and our expectations for scaling near 
a saddle-node bifurcation. 

We studied the effects of the rewiring on the network 
geometry. Degree distributions were altered, and mean 
field arguments were able to predict the distributions for 
susceptibles and recovereds. However, a new analytical 
approach that includes correlations is needed to fully un- 
derstand the degree distributions. 

The other network property we considered was the 
distribution of distances from non-infected nodes to the 
nearest infective, a quantity that may be important in 
disease spreading and its control. This distribution de- 
pended primarily on node and link dynamics rather than 
on higher order correlations, so it could be predicted from 
random graphs. It is possible to compute the distribu- 
tion of distances from an infective analytically for ran- 
dom graphs, but this calculation is awkward for a three 
species system (S, I, R) and does not have a simple func- 
tional form, so we have omitted discussion of analytical 
results here. 

Delayed outbreaks were observed in the rewired sys- 
tem. Peaks in the infective fraction lagged behind peaks 
in the number of nodes that neighbor an infective. For 
the parameter values studied here, the lag time is on the 
order of 10% of the mean infectious period, which might 
be considered a very short lag time. However, the pa- 
rameters in this study were not selected to correspond 
to any specific disease. For most real diseases, we would 
expect a much slower rate for immunity to wear off and 
recovered individuals to become susceptible again, com- 
pared to the mean infectious period. This regime would 
be more difficult to study in Monte Carlo simulations, 
since the average number of infected nodes would be 
much smaller than seen here. Further work is needed to 
determine whether the observation of delayed outbreaks 
due to rewiring would persist or perhaps become more 
significant in a physically realistic system. 

Finally, we considered lifetimes of the endemic state 
near the saddle-node bifurcation where it loses stability. 
In order to achieve extinction from a steady state, the dis- 
ease must first overcome the attractive forces, which are 
weak near the bifurcation point. Due to the generic local 
topology of the saddle-node structure, the escape rate is 
well characterized analytically j2^]. We found that the 



Monte Carlo simulation agrees qualitatively with the es- 
cape times and yields a well-known power law. However, 
this is for parameters in which endemic and extinction 
states are not too far apart. Such extinction regimes 
can be analyzed using a Fokker-Planck approach [13]. 
Although preliminary results are in agreement with the 
expected scaling, further study of the lifetime scaling is 
needed, including in regimes that are physically realistic, 
and where the usual extinction rates cannot be modeled 
with a Fokker-Planck approach. 

In addition to the directions for further research men- 
tioned above, a major challenge is to develop network 
geometries and rewiring rules that are more consistent 
with human social networks. Real social networks are 
expected to have community structure (28| . The net- 
works studied here did not. In our rewired networks, 
the number of connections from non-infected nodes to 
infected nodes was reduced in comparison to a random 
network, while the number of connections between non- 
infected nodes was increased. However, this process does 
not induce a community structure on the network in the 
Newman-Girvan Q-modularity sense [29| . Rewiring was 
a non-local effect; new neighbors were chosen at random 
amongst all non-infected nodes in the network rather 
than introducing a local structure. It has been shown for 
static networks that community structure affects disease 
dynamics (sol . [3l| , and we expect an impact in adaptive 
networks as well. 

In the current work, an ideal setting was proposed 
where non-infected nodes were assumed to behave ra- 
tionally and have perfect knowledge of the disease status 
of their current neighbors and potential new neighbors. 
It would be of interest to consider a situation in which 
not all contagious individuals appear ill or know that 
they arc contagious, as might be the case for a sexu- 
ally transmitted disease possessing asymptomatic indi- 
viduals.. This effect might be modeled by simultaneous 
spreading through the network of both the disease and 
information about the disease. If the current model is ex- 
tended with information and community structure, social 
dynamics could be extrapolated to improve contact trac- 
ing and epidemic control in organized populations with 
local structure. 
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